テーブルデータの特異値分解

はじめに

　テーブルデータを行列とみなし、各種の線形代数手法を適用することで、データの特徴の可視化、潜在因子の抽出、容量圧縮などができる。その結果は分析者個人に依らない。ならば1変数、2変数のデータから基本的な統計量、分布、散布図を示すことがマナーであるように、多変数のテーブルデータから、その相関行列、主成分分析または特異値分解を示すことも基本的な記述統計ではないかと考える。データ分析の前段階としてそれらの十分な記述統計を準備することで、データから仮説検証するときも、予測モデルをつくるときも、（経験、勘でなく）その分析方法をより科学的な根拠から説明することができるだろう。

　ここではテーブルデータの特異値分解の基本についてまとめる。特異値分解は画像分析、音声分析、リコメンデーションなど実社会のビッグデータに対して最も適用されている分析手法である。今回は小さなデータを用いた説明になるが、小さな行列で成立する方法をそのまま何万次元のデータにも適用できることが線形代数手法の長所でもある。なお特異値分解を含め、行列分解の概要については別の小文でもまとめている。

(https://kng-mtd.github.io/test0/matrixDecomposition)

特異値分解、主成分分析は直交する因子を抽出する

grViz("digraph {
  graph [label='観測値と因子', labelloc=t,
    layout = dot, rankdir = TB,
    nodesep=0.5,ranksep=1]

  
  node [shape=rectangle]
  f1 [label = '体力',shape=ellipse]
  f2 [label = '知力',shape=ellipse]
  f3 [label = '運',shape=ellipse]
  v1 [label = 'マージャン']
  v2 [label = '年俸']  
  v3 [label = '100m走']
  v4 [label = '遠投']

  f1->v1
  f1->v2
  f1->v3
  f1->v4
  f2->v1
  f2->v2
  f2->v3
  f2->v4
  f3->v1
  f3->v2
  f3->v3
  f3->v4

}")

選手	100ｍ走	遠投	マージャン	年棒
しおみ ながおか やまだ むらかみ オスナ サンタナ なかむら やまざき いしかわ	6.0 6.1 6.5 7.0 7.7 7.5 7.5 6.5 7.7	100 80 70 70 60 50 80 100 70	3 4 6 4 8 6 8 4 3	10,000 5,000 20,000 20,000 15,000 15,000 12,000 3,000 7,000

　たとえば、スポーツ選手を比較しようとするとき、まずいろいろな観測値を集める。それである程度の各選手の特徴を把握できる。しかし観測値は種類が多く（いくらでも集められる）、また観測値は多様な因子の組み合わせで決まるので、ある観測値が大きいと別の観測値も大きい、小さいなどの複雑な関係がある。分析者により、どの観測値を選び、どのように組み合わせるかに差、クセがでてくる、「いや、オレはこの数字とこの数字を見るよ」。

　もっと一般的に、体系的に対象の特徴抽出、比較、分類するには、複雑な関係にある観測値ではなく、観測値の元となる、かつ互いに関係を持たない因子を見つけるといい。因子間に関係がなければ、単純に個々の因子を大きい順に並べて、わかりやすく比較することができる。各因子が互いに関係を持たないとき（無相関のとき）、それらを「互いに直交する」という。主成分分析、特異値分解はそのような直交する因子をみつけて、対象を比較、分類、可視化する手法である（データ特有のハイパーパラメータなども不要）。

　特異値分解は、データ行列から互いに直交する因子を抽出し、各因子の重要性、各対象の特徴、各変数の特徴を同時に可視化するものである。対象を行、その変数を列とするデータ行列を、対象の特徴を表す直交する因子（左特異ベクトル）、変数の特徴を表す直交する因子（右特異ベクトル）、因子の重要性を表す特異値行列に分解するもの。

　主成分分析は、データがもつ変数の相関行列から互いに直交する因子（主成分）を抽出し、各因子の重要性、各対象の因子の大きさ（主成分スコア）を可視化する。 　 　 　変数間の相関行列を固有値分解し、固有ベクトル＝主成分と、固有値＝主成分の重要性、各対象がもつ主成分方向の大きさ（主成分スコア）を計算する。数学的には、互いに直行する分散の大きい方向をみつけることは、相関行列の固有値分解問題であり、対称行列である相関行列がもつ固有ベクトルは直交することが示されている。

　前処理により特異値分解、主成分分析は対象、変数の結果は同じになる（どちらでもいい）。計算が早い特異値分解、結果を理解しやすい主成分分析あたりで使い分けられている。

データ行列の特異値分解

データとその前処理

　例として３つの島に生息する、３種のペンギンの性別、体重、羽の長さ、嘴の長さ、嘴の幅の3年間調査であるpalmerpenguinsデータを使う。

tb0=penguins
glimpse(tb0)

## Rows: 344
## Columns: 8
## $ species           <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island            <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm     <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g       <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex               <fct> male, female, female, NA, female, male, female, male…
## $ year              <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…

summary(tb0)

##       species          island    bill_length_mm bill_depth_mm 
##  Adelie   :152   Biscoe   :168   Min.   :32.1   Min.   :13.1  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.2   1st Qu.:15.6  
##  Gentoo   :124   Torgersen: 52   Median :44.5   Median :17.3  
##                                  Mean   :43.9   Mean   :17.1  
##                                  3rd Qu.:48.5   3rd Qu.:18.7  
##                                  Max.   :59.6   Max.   :21.5  
##                                  NA's   :2      NA's   :2     
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172       Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190       1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197       Median :4050   NA's  : 11   Median :2008  
##  Mean   :201       Mean   :4202                Mean   :2008  
##  3rd Qu.:213       3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231       Max.   :6300                Max.   :2009  
##  NA's   :2         NA's   :2

#datatable(tb0,options=list(scrollX=500))

欠損値のある対象を除く。年度を量的変数から質的変数に型変換する。

tb0 %<>%
  drop_na() %>% 
  mutate(year=as.factor(year))

量的変数の正規化

　特異値分解でも主成分分析と同様、変数の[0,1]正規化や標準正規化をする、変数の単位が異なるときは必須。一部の変数を正規化することが障害になるかどうか、特異値分解前の変数の記述統計により検討し、その変数を除くことで対処する。もし正規化しない場合、特異値分解の結果の解釈を難しくする。

tbn=tb0 %>%
  select(is.numeric) # just numeric variable

nn=names(tbn)
l=c()
u=c()

for(i in nn){
  l[i]=min(tbn[i])
  u[i]=max(tbn[i])
  tbn[i]=scale(tbn[i],center=l[i],scale=u[i]-l[i])
}

質的変数のダミー変数化

　質的変数を含むデータに対しても線形代数手法を適用するために、２値の質的変数は{0,1}のダミー変数、3水準以上の質的変数は複数の{0,1}のダミー変数に変換する。また質的変数が順序をもつとき、整数値の連続変数に変換することでほぼ対処できる。

factor2ind=function(x, baseline){
  xname=deparse(substitute(x))
  n=length(x)
  x=as.factor(x)
  if(!missing(baseline)) x=relevel(x, baseline)
  X=matrix(0L, n, length(levels(x)))
  X[(1:n) + n*(unclass(x)-1)]=1L
  X[is.na(x),]=NA
  dimnames(X)=list(names(x), paste(xname, levels(x), sep = ":"))
  return(X[,-1,drop=FALSE])
}

tbc0=tb0 %>%
  select(!is.numeric) # just categorical variable

nc=names(tbc0)

tbc=tibble(rep(NA,nrow(tbc0)))
for(i in nc){
  btb=as_tibble(factor2ind(tbc0[[i]]))
  names(btb)=names(btb) %>%
    str_replace('tbc0\\[\\[i\\]\\]',i)
  tbc=tbc %>%
    bind_cols(btb)
}
tbc=tbc[,-1]

　すべての変数から成る相関行列をみて、他の変数との相関が微小な連続変数を除く（当然、質的変数の別水準を表すダミー変数と逆相関があるだけのダミー変数も除く）。そのような変数は特異値分解の結果に影響しない。

tb=bind_cols(tbn,tbc)

par(mfrow=c(1,1))
corrplot(cor(tb))

tb=tb %>% 
  select(-'year:2008',-'year:2009')
tbc=tbc %>% 
  select(-'year:2008',-'year:2009')

特異値分解する、前処理後のデータ行列

glimpse(tb)

## Rows: 333
## Columns: 9
## $ bill_length_mm      <dbl> 0.2545, 0.2691, 0.2982, 0.1673, 0.2618, 0.2473, 0.…
## $ bill_depth_mm       <dbl> 0.667, 0.512, 0.583, 0.738, 0.893, 0.560, 0.774, 0…
## $ flipper_length_mm   <dbl> 0.1525, 0.2373, 0.3898, 0.3559, 0.3051, 0.1525, 0.…
## $ body_mass_g         <dbl> 0.292, 0.306, 0.153, 0.208, 0.264, 0.257, 0.549, 0…
## $ `species:Chinstrap` <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `species:Gentoo`    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `island:Dream`      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `island:Torgersen`  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,…
## $ `sex:male`          <int> 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1,…

summary(tb)

##  bill_length_mm  bill_depth_mm   flipper_length_mm  body_mass_g   
##  Min.   :0.000   Min.   :0.000   Min.   :0.000     Min.   :0.000  
##  1st Qu.:0.269   1st Qu.:0.298   1st Qu.:0.305     1st Qu.:0.236  
##  Median :0.451   Median :0.500   Median :0.424     Median :0.375  
##  Mean   :0.432   Mean   :0.484   Mean   :0.491     Mean   :0.419  
##  3rd Qu.:0.600   3rd Qu.:0.667   3rd Qu.:0.695     3rd Qu.:0.576  
##  Max.   :1.000   Max.   :1.000   Max.   :1.000     Max.   :1.000  
##  species:Chinstrap species:Gentoo   island:Dream   island:Torgersen
##  Min.   :0.000     Min.   :0.000   Min.   :0.000   Min.   :0.000   
##  1st Qu.:0.000     1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.000   
##  Median :0.000     Median :0.000   Median :0.000   Median :0.000   
##  Mean   :0.204     Mean   :0.357   Mean   :0.369   Mean   :0.141   
##  3rd Qu.:0.000     3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:0.000   
##  Max.   :1.000     Max.   :1.000   Max.   :1.000   Max.   :1.000   
##     sex:male    
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :1.000  
##  Mean   :0.505  
##  3rd Qu.:1.000  
##  Max.   :1.000

特異値分解の結果

　特異値分解では、データ行列（行が各対象、列は変数）の元の変数と同じ数（本データは９コ）の特異値が計算される。 適当な「次数」＝因子数を決めて、特異値分解を打ち切る。対象数×次数の左特異ベクトル行列、次数次元の正方行列である特異値行列、変数の数×次数の右特異ベクトル行列ができる。

　たとえば次数を４としたとき、左特異ベクトル行列は対象を４コの数値（＝因子の大きさ）で表したもの。特異値行列は各因子が対象の特徴をどれだけ表すかの重み。右特異ベクトル行列は各因子が元の各変数とどれだけ相関するかを示す重みであり、すなわち元の変数の特徴を４コの数値で表したものとなる。

k=4
ku=k
kv=k
s=svd(tb,ku,kv) #u:n x ku, d:ku x kv, v:m x kv
#s
#s$d #singular value

u=s$u #left singular matrix
head(u,10)

##          [,1]    [,2]   [,3]    [,4]
##  [1,] -0.0458 0.00767 0.1258 -0.0539
##  [2,] -0.0243 0.00181 0.0571 -0.1394
##  [3,] -0.0262 0.00308 0.0569 -0.1454
##  [4,] -0.0269 0.00615 0.0638 -0.1493
##  [5,] -0.0517 0.01009 0.1297 -0.0669
##  [6,] -0.0225 0.00498 0.0603 -0.1386
##  [7,] -0.0555 0.00179 0.1239 -0.0661
##  [8,] -0.0219 0.00592 0.0585 -0.1383
##  [9,] -0.0533 0.01063 0.1316 -0.0701
## [10,] -0.0554 0.00613 0.1314 -0.0717

d=diag(s$d[1:k])
d

##      [,1] [,2] [,3] [,4]
## [1,] 23.1  0.0 0.00 0.00
## [2,]  0.0 13.1 0.00 0.00
## [3,]  0.0  0.0 8.94 0.00
## [4,]  0.0  0.0 0.00 6.74

v=s$v #right singular matrix
head(v,10)

##          [,1]    [,2]    [,3]    [,4]
##  [1,] -0.3620 -0.0313 -0.1654 -0.0997
##  [2,] -0.3633  0.2616  0.2438 -0.2756
##  [3,] -0.4026 -0.2163 -0.1358 -0.1665
##  [4,] -0.3501 -0.2091 -0.0413 -0.0497
##  [5,] -0.1816  0.3939 -0.3387 -0.0689
##  [6,] -0.3078 -0.5738 -0.3605 -0.0478
##  [7,] -0.2911  0.5961 -0.3272  0.0387
##  [8,] -0.0762  0.0134  0.4749 -0.7163
##  [9,] -0.4838  0.0147  0.5620  0.6019

　元の行列は左特異ベクトル行列、特異値行列、右特異ベクトル行列の積で近似される（近似行列）。特異値分解による近似行列は、あらゆる行列分解のうち、元の行列に最も近似する（行列要素間の差のﾌﾛﾍﾞﾆｳｽﾉﾙﾑ最小）（Eckart-Youngの定理）。左右の特異ベクトル行列、特異値行列の要素数の合計が、元のデータ行列の要素数よりも少なくなれば圧縮できるが、どれだけ元の行列と近似行列が近いのかを確認するには、以下のような方法があるだろう。

• 行列全体の数値をヒートマップで（直感的に）比べる

• 量的変数の各要素の差をみる（平均二乗誤差 MSE、二乗誤差の分散 V[SE]）

• 近似行列の各量的変数の基本統計量をみる

• 各質的変数のクロス表をみる（近似行列の数値は0.5を境界に0,1に分ける）

• 各量的変数の分布をボックスプロット、ヒストグラムで（直感的に）比べる

• 各量的変数の平均の差をｔ検定する

• 相関行列を比べる

• 散布図行列を（直感的に）比べる

row.names(v)=names(tb)

tb_svd0=as_tibble(u %*% d %*% t(v))
tbn_svd=tb_svd0 %>%
  select(names(tbn))
tbc_svd=tb_svd0 %>%
  select(names(tbc))
tbc_svd[tbc_svd<.5]=0
tbc_svd[tbc_svd>.5]=1

tb_svd=bind_cols(tbn_svd,tbc_svd)

heatmap(as.matrix(tb),Rowv=NA,Colv=NA,scale='none',margins=c(10,3),main='origin')

#heatmap(u %*% d %*% t(v),Rowv=NA,Colv=NA,scale='none',margins=c(10,3))
heatmap(as.matrix(tb_svd),Rowv=NA,Colv=NA,scale='none',margins=c(10,3),main='SVD4')

cat('MSE: ',mean((tbn-tbn_svd)**2 |> unlist()))

## MSE:  0.00859

cat('V[SE]: ',var((tbn-tbn_svd)**2 |> unlist()))

## V[SE]:  0.000335

summary(tb_svd)

##  bill_length_mm  bill_depth_mm   flipper_length_mm  body_mass_g   
##  Min.   :0.108   Min.   :0.039   Min.   :0.107     Min.   :0.051  
##  1st Qu.:0.303   1st Qu.:0.301   1st Qu.:0.292     1st Qu.:0.217  
##  Median :0.478   Median :0.452   Median :0.417     Median :0.383  
##  Mean   :0.439   Mean   :0.462   Mean   :0.489     Mean   :0.415  
##  3rd Qu.:0.576   3rd Qu.:0.627   3rd Qu.:0.705     3rd Qu.:0.563  
##  Max.   :0.759   Max.   :0.957   Max.   :0.996     Max.   :0.894  
##  species:Chinstrap species:Gentoo   island:Dream   island:Torgersen
##  Min.   :0.000     Min.   :0.000   Min.   :0.000   Min.   :0.000   
##  1st Qu.:0.000     1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.000   
##  Median :0.000     Median :0.000   Median :0.000   Median :0.000   
##  Mean   :0.216     Mean   :0.357   Mean   :0.369   Mean   :0.141   
##  3rd Qu.:0.000     3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:0.000   
##  Max.   :1.000     Max.   :1.000   Max.   :1.000   Max.   :1.000   
##     sex:male    
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :1.000  
##  Mean   :0.505  
##  3rd Qu.:1.000  
##  Max.   :1.000

for(i in names(tbc)){
  cat('\n',i)
  table(tbc[[i]],tbc_svd[[i]]) %>% print()
  #wilcox.test(tb[[i]],tb_svd[[i]]) %>% print()
}

## 
##  species:Chinstrap   
##       0   1
##   0 261   4
##   1   0  68
## 
##  species:Gentoo   
##       0   1
##   0 214   0
##   1   0 119
## 
##  island:Dream   
##       0   1
##   0 210   0
##   1   0 123
## 
##  island:Torgersen   
##       0   1
##   0 286   0
##   1   0  47
## 
##  sex:male   
##       0   1
##   0 165   0
##   1   0 168

par(mfrow=c(1,1))
boxplot(tbn,main='origin')

boxplot(tbn_svd,main='SVD4')

par(mfrow=c(2,1))
for(i in names(tbn)){
  hist(tb[[i]],main=i,xlim=c(0,1.5))
  hist(tb_svd[[i]],main='',xlim=c(0,1.5))
  
  #a=cor(tb[i],tb_svd[i])
  #cat(i,' correlation:',a[[1]],'\n')
  t.test(tbn[i],tbn_svd[i]) %>% print()
}

## 
##  Welch Two Sample t-test
## 
## data:  tbn[i] and tbn_svd[i]
## t = -0.4, df = 638, p-value = 0.7
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0339  0.0213
## sample estimates:
## mean of x mean of y 
##     0.432     0.439

## 
##  Welch Two Sample t-test
## 
## data:  tbn[i] and tbn_svd[i]
## t = 1, df = 663, p-value = 0.2
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0140  0.0587
## sample estimates:
## mean of x mean of y 
##     0.484     0.462

## 
##  Welch Two Sample t-test
## 
## data:  tbn[i] and tbn_svd[i]
## t = 0.1, df = 664, p-value = 0.9
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0334  0.0381
## sample estimates:
## mean of x mean of y 
##     0.491     0.489

## 
##  Welch Two Sample t-test
## 
## data:  tbn[i] and tbn_svd[i]
## t = 0.2, df = 664, p-value = 0.8
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0301  0.0373
## sample estimates:
## mean of x mean of y 
##     0.419     0.415

cat('origin')

## origin

cor(tbn)

##                   bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm             1.000        -0.229             0.653       0.589
## bill_depth_mm             -0.229         1.000            -0.578      -0.472
## flipper_length_mm          0.653        -0.578             1.000       0.873
## body_mass_g                0.589        -0.472             0.873       1.000

cat('SVD4')

## SVD4

cor(tbn_svd)

##                   bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm            1.0000       -0.0389             0.817       0.763
## bill_depth_mm            -0.0389        1.0000            -0.383      -0.304
## flipper_length_mm         0.8172       -0.3826             1.000       0.972
## body_mass_g               0.7630       -0.3038             0.972       1.000

plot(tbn,main='origin')

plot(tbn_svd,main='SVD4')

# psyck::r.test
# 1.H0 r12=0, n single group
#   r.test(n,r12)
# 2,H0 r12=r34, n1,n2 not correspond groups
#   r,test(n=n1,r12,n2=n2,r34)
# 3,H0 r12=r13, n correspond groups
#   r.test(n,r12,r23,r13)
# 4.H0 r12=r34, n correspond groups
#   r.test(n,r12,r34,r23,r13,r14,r24)

#rtest=r.test(nrow(tbn),cor(tbn),cor(tbn_svd))
#rtest$p

　近似行列において、量的変数の平均、中央値は再現されている。ダミー変数の値{0,1}を、近似行列では{0,1}に近い数値でほぼ再現されている。分布の形はほぼ再現されるが、ノイズ、外れ値が除去されている。相関関係の大小がより強調されている。散布図から、複数のクラスの別々の相関関係をみることができる。

特異値分解の結果から新たな「観測値」を予測する

　観測値を単純化した近似行列を、母集団の母数とすることができるかもしれない。観測値はそれにノイズが加わったものとするモデルが考えられる。そこで特異値分解の結果にノイズを加えた「観測値」をつくってみる。

　近似行列の全要素に上下10%の一様乱数を加えたものを10コから成る、対象数を10倍にした拡張データをつくった。拡張データから元と同じ対象数をサンプリングをすることで「観測」をシミュレーションした。

expd0=tibble()
for(i in 1:10){ # 10times
  tb1=as_tibble(u %*% d %*% t(v))*runif(nrow(tb)*ncol(tb),0.9,1.1) # ±10%
  expd0=bind_rows(expd0,tb1)
}

expdn=expd0 %>%
  select(names(tbn))
expdc=expd0 %>%
  select(names(tbc))
expdc[expdc<.5]=0
expdc[expdc>.5]=1

expd=bind_cols(expdn,expdc)

#glimpse(expd)
datatable(expd,options=list(scrollX=500))

summary(expd)

##  bill_length_mm  bill_depth_mm   flipper_length_mm  body_mass_g   
##  Min.   :0.099   Min.   :0.037   Min.   :0.102     Min.   :0.047  
##  1st Qu.:0.296   1st Qu.:0.295   1st Qu.:0.291     1st Qu.:0.217  
##  Median :0.470   Median :0.452   Median :0.415     Median :0.381  
##  Mean   :0.439   Mean   :0.462   Mean   :0.488     Mean   :0.415  
##  3rd Qu.:0.572   3rd Qu.:0.630   3rd Qu.:0.709     3rd Qu.:0.573  
##  Max.   :0.820   Max.   :1.043   Max.   :1.060     Max.   :0.966  
##  species:Chinstrap species:Gentoo   island:Dream   island:Torgersen
##  Min.   :0.000     Min.   :0.000   Min.   :0.000   Min.   :0.000   
##  1st Qu.:0.000     1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.000   
##  Median :0.000     Median :0.000   Median :0.000   Median :0.000   
##  Mean   :0.231     Mean   :0.357   Mean   :0.369   Mean   :0.141   
##  3rd Qu.:0.000     3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:0.000   
##  Max.   :1.000     Max.   :1.000   Max.   :1.000   Max.   :1.000   
##     sex:male    
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :1.000  
##  Mean   :0.505  
##  3rd Qu.:1.000  
##  Max.   :1.000

write_csv(expd,'penguins_expand.csv')

expd=read_csv('penguins_expand.csv')

tb1=sample_n(expd,333) #sample n=333 
tb1n=tb1 %>% select(names(tbn))

par(mfrow=c(1,1))
boxplot(tbn,main='origin')

boxplot(tb1n,main='new obs.')

par(mfrow=c(2,1))
for(i in names(tbn)){
  hist(tb[[i]],main=i,xlim=c(0,1.5))
  hist(tb1[[i]],main='',xlim=c(0,1.5))

  t.test(tb[i],tb1[i]) %>% print()
}

## 
##  Welch Two Sample t-test
## 
## data:  tb[i] and tb1[i]
## t = -0.5, df = 636, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0345  0.0206
## sample estimates:
## mean of x mean of y 
##     0.432     0.439

## 
##  Welch Two Sample t-test
## 
## data:  tb[i] and tb1[i]
## t = 0.6, df = 658, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0264  0.0486
## sample estimates:
## mean of x mean of y 
##     0.484     0.473

## 
##  Welch Two Sample t-test
## 
## data:  tb[i] and tb1[i]
## t = 0.8, df = 661, p-value = 0.4
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02  0.05
## sample estimates:
## mean of x mean of y 
##     0.491     0.476

## 
##  Welch Two Sample t-test
## 
## data:  tb[i] and tb1[i]
## t = 0.6, df = 660, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0233  0.0424
## sample estimates:
## mean of x mean of y 
##     0.419     0.409

cat('origin')

## origin

cor(tbn)

##                   bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm             1.000        -0.229             0.653       0.589
## bill_depth_mm             -0.229         1.000            -0.578      -0.472
## flipper_length_mm          0.653        -0.578             1.000       0.873
## body_mass_g                0.589        -0.472             0.873       1.000

cat('new obs.')

## new obs.

cor(tb1n)

##                   bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm            1.0000        0.0349             0.785       0.738
## bill_depth_mm             0.0349        1.0000            -0.353      -0.277
## flipper_length_mm         0.7848       -0.3526             1.000       0.951
## body_mass_g               0.7379       -0.2773             0.951       1.000

plot(tbn,main='origin')

plot(tb1n,main='new obs.')

ランダム行列による特異値分解

library(imager)

img0=load.image('image1.jpg')
#plot(img0,axes=F)

img=grayscale(img0) #[0,1]
#plot(img,axes=F)
save.image(img,'image10.jpg')

mx=as.matrix(img)
img=as.cimg(mx)
plot(img,axes=F,main='origin')

w0=dim(mx)[1]
h0=dim(mx)[2]
cat('width: ',w0,'   height: ',h0)

## width:  5472    height:  3648

img1=resize(img,w0/64,h0/64) 
plot(img1,axes=F,main='1/64')

k=20

gc(reset=T)

##             used (Mb) gc trigger (Mb)  max used (Mb)
## Ncells   2554552  136    4124740  220   2554552  136
## Vcells 144330088 1101  267872290 2044 144330088 1101

#cat('\nSVD k:',k)
system.time({
  s=svd(mx,k,k)
  u=s$u
  d=diag(s$d[1:k])
  v=s$v
  tb1=u %*% d %*% t(v)
})  

##   ユーザ システム     経過 
##    92.51     0.64   195.02

  img=as.cimg(tb1)
  plot(img,axes=F,main=str_c('SVD ',k))

  save.image(img,'image1svd20.jpg')

  
gc(reset=T)

##             used (Mb) gc trigger (Mb)  max used (Mb)
## Ncells   2555711  136    4124740  220   2555711  136
## Vcells 184443856 1407  329226953 2512 184443856 1407

#cat('\nrSVD k:',k)
system.time({
  s=rsvd(mx,k)
  u=s$u
  d=diag(s$d[1:k])
  v=s$v
  tb1=u %*% d %*% t(v)
})

##   ユーザ システム     経過 
##     1.39     0.02     3.41

  img=as.cimg(tb1)
  plot(img,axes=F,main=str_c('rSVD ',k))

  save.image(img,'image1rsvd20.jpg')



k=50

gc(reset=T)

##             used (Mb) gc trigger (Mb)  max used (Mb)
## Ncells   2557223  137    4124740  220   2557223  137
## Vcells 224366754 1712  395152343 3015 224366754 1712

#cat('\nSVD k:',k)
system.time({
  s=svd(mx,k,k)
  u=s$u
  d=diag(s$d[1:k])
  v=s$v
  tb1=u %*% d %*% t(v)
})

##   ユーザ システム     経過 
##    72.92     0.38   184.48

  img=as.cimg(tb1)
  plot(img,axes=F,main=str_c('SVD ',k))

  save.image(img,'image1svd50.jpg')

  
gc(reset=T)

##             used (Mb) gc trigger (Mb)  max used (Mb)
## Ncells   2557437  137    4124740  220   2557437  137
## Vcells 244612716 1866  404516612 3086 244612716 1866

#cat('\nrSVD k:',k)
system.time({
  s=rsvd(mx,k)
  u=s$u
  d=diag(s$d[1:k])
  v=s$v
  tb1=u %*% d %*% t(v)
})

##   ユーザ システム     経過 
##     2.52     0.10     6.03

  img=as.cimg(tb1)
  plot(img,axes=F,main=str_c('rSVD ',k))

  save.image(img,'image1rsvd50.jpg')